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Abstract 

Deformation  twinning  at  the  tip  of  a  straight  crack  or  notch  is  analyzed  using  a  phase-field  method  that  seeks  equilibrium  twin  mor¬ 
phologies  via  direct  minimization  of  a  free  energy  functional.  For  isotropic  solids,  the  tendency  to  twin  under  mode  I  or  mode  II  loading 
is  found  to  depend  weakly  on  Poisson’s  ratio  and  elastic  nonlinearity  and  strongly  on  surface  energy  and  twinning  shear  (i.e.  eigenstrain). 
Depending  on  the  coherent  twin  boundary  energy,  anisotropy  of  surface  energy  is  important  for  mode  I  loading  but  less  so  for  mode  II. 
Model  predictions  for  several  single  crystals  are  in  agreement  with  experimental  observations.  Calcite  demonstrates  a  preference  for 
mode  I  cleavage  crack  extension  over  crack  tip  twinning.  Magnesium  shows  a  likelihood  for  tensile  twinning  from  a  pre-existing  crack 
on  the  basal  plane.  In  sapphire,  a  preference  for  rhombohedral  twins  over  basal  twins  is  apparent,  with  the  latter  thinner  in  shape  than 
the  former. 

Published  by  Elsevier  Ltd.  on  behalf  of  Acta  Materialia  Inc. 

Keywords:  Phase  field;  Twinning;  Fracture;  Crystals;  Modeling  and  simulation 


1.  Introduction 

Deformation  twinning,  i.e.  twinning  induced  by 
mechanical  stress,  and  cleavage  fracture,  i.e.  transgranular 
fracture  on  preferred  crystallographic  planes,  are  two  fun¬ 
damental  inelastic  deformation  mechanisms  that  may 
occur  in  crystals.  In  some  cases,  twins  or  twin  boundaries 
may  act  as  nucleation  sites  for  fracture  [1];  in  others,  crack 
tips  may  provide  the  necessary  stress  concentrations  for 
twin  nucleation  and  growth  [2]. 

In  dynamic  experiments  such  as  plate  impact,  both  twin¬ 
ning  and  cleavage  fracture  may  occur,  but  whether  twin¬ 
ning  precedes  or  follows  fracture  cannot  usually  be 
determined  from  analysis  of  experimental  (e.g.  Hugoniot) 
data  and  post-mortem  material  characterization.  This  is 
the  case  for  impact  experiments  on  sapphire  [3],  wherein 
theoretical  strengths  for  twinning,  slip  and  shear  fracture 
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are  thought  to  be  of  the  same  order  of  magnitude  [4]. 
Indentation  experiments  on  ceramics  and  minerals  often 
show  evidence  of  both  fracture  and  twinning,  with  surface 
damage  strongly  promoting  the  nucleation  and  growth  of 
twins  in  calcite,  for  example  [5,6]. 

The  present  work  seeks  (i)  to  develop  an  improved 
understanding  of  the  general  factors  affecting  twinning 
induced  by  fracture;  and  (ii)  to  test  a  predictive  model  for 
twinning  at  a  crack/notch  tip  in  several  real  materials. 
Phase-field  theory  and  numerical  simulation  are  applied 
to  study  twin  nucleation  and  growth  from  a  pre-existing 
crack  or  notch.  A  prior  analysis  [2]  based  on  the  Peierls- 
Nabarro  concept  and  ideas  from  Ref.  [7]  was  developed 
to  judge  the  tendency  of  a  solid  to  undergo  either  microt- 
winning  or  slip  of  leading  and  trailing  partial  dislocations 
on  the  same  plane.  This  treatment  demonstrated  success 
for  several  face-centered  cubic  metals  when  compared  to 
results  of  atomic  simulations,  but  the  model  requires 
knowledge  of  parameters  associated  with  the  stacking  fault 
energy  surface  that  seem  only  to  be  obtainable  from  atomic 
simulation  of  planar  defects.  An  analytical  model 
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predicting  the  likelihood  of  twinning  or  extension  of  a 
mode  I  slit  crack  is  described  in  Ref.  [8];  this  model  requires 
knowledge  of  energetic  data  associated  with  total  resistance 
to  twinning  or  fracture  that  are  evidently  not  available 
from  standard  experiments,  though  twin  boundary  and 
cleavage  surface  energies  can  be  substituted  as  an  approx¬ 
imation.  As  will  be  shown  later,  the  current  work  offers 
more  insight  into  the  total  twinning  resistance  that  could 
be  used  in  such  a  model. 

In  the  present  phase-field  approach  [9,10],  the  only  con¬ 
stitutive  model  parameters  are  the  twinning  shear  (usually 
known  from  crystallography),  elastic  constants  (known 
from  experiments),  gradient  energy  coefficient(s),  and  dou¬ 
ble-well  energy  barrier  height.  The  latter  two  can  be  related 
to  the  twin  boundary  surface  energy  (measurable  from 
experiment  or  atomic  simulation)  and  the  thickness  of  the 
diffuse  interface,  which  can  be  assigned  physical  signifi¬ 
cance  if  the  normal  distance  from  the  phase/twin  boundary 
over  which  atoms  are  displaced  from  their  usual  positions 
in  a  perfect  crystal  is  known,  e.g.  from  atomic  simulation. 
For  prescribed  boundary  conditions  and  problem  geome¬ 
try,  the  tendency  for  a  pre-cracked  crystal  to  continue  to 
crack  or  to  twin  then  necessarily  depends  on  these  param¬ 
eters  and  the  surface  energy  associated  with  cleavage  frac¬ 
ture.  The  model  does  not  depend  on  the  numerical  method 
of  solution  or  grid  size  (i.e.  the  theory  itself  is  mesh  inde¬ 
pendent),  but  any  discretization  should  be  sufficiently 
refined  to  resolve  continuum  fields  where  gradients  exist, 
e.g.  in  interfacial  regions. 

The  present  work  does  not  address  plastic  slip  distinct 
from  the  motion  of  twinning  partials  inherent  in  deforma¬ 
tion  twinning.  Nanoscale  treatments  of  discrete  slip  criteria 
include  Refs.  [2,7];  mesoscale  continuum  crystal  plasticity 
models  of  slip  and  twinning  are  described  in  Ref.  [11] 
and  references  therein. 

Several  other  relevant  modeling  approaches  are  noted. 
Phase  transformation  has  been  studied  at  tips  of  moving 
cracks  via  analytical  solutions  to  phase-field  models 
[12,13].  Like  twinning,  phase  transformations  may  be 
induced  by  strong  local  elastic  fields  at  crack  tips;  transfor¬ 
mation  strain  in  the  former  often  includes  dilatation,  while 
deformation  twinning  involves  large  shear  without  inelastic 
volume  change.  Competition  among  phase  transformation, 
fracture  and  plastic  deformation  was  studied  using  a  con¬ 
tinuum  thermodynamic  approach  implemented  in  the 
finite-element  method  [14];  twinning  was  also  modeled. 
Phase-field  models  of  fracture  have  also  been  implemented 
[15-19].  The  present  paper  does  not  develop  a  phase-field 
model  for  fracture — herein  a  stationary  pre-crack  is  repre¬ 
sented  explicitly  by  a  thin  notch  with  free  surface  boundary 
conditions — but  conceivably  both  twinning  and  fracture 
could  be  modeled  simultaneously  using  the  phase-field 
approach,  with  distinct  order  parameters  accounting  for 
transformation  to  twinned  and/or  fractured  material. 

This  paper  is  organized  as  follows.  The  phase-field  the¬ 
ory  developed  and  implemented  in  Refs.  [9,10]  is  reviewed 
in  Section  2,  including  various  elasticity  models  (linear 


isotropic,  linear  anisotropic,  nonlinear  neo-Hookean)  con¬ 
sidered  later.  Analysis  of  possible  twinning  or  crack  exten¬ 
sion  under  pure  mode  I  or  mode  II  loading  in  generic 
isotropic  elastic  solids  follows  in  Section  3.  In  this  analysis, 
a  normalized  energy  functional  is  derived  that  depends  on 
several  dimensionless  material  parameters.  The  effects  of 
these  parameters — twinning  shear,  Poisson’s  ratio  and 
normalized  twin  boundary  energy — on  twinnability  are 
investigated  through  phase-field  simulations.  In  Section  4, 
two-dimensional  (2-D)  simulations  of  twinning  from  a 
mode  I  crack  are  reported  for  calcite  (CaC03),  sapphire 
(a-Al203)  and  magnesium  (Mg)  and  compared  with  exper¬ 
imental  observations.  In  Section  5,  three-dimensional  (3-D) 
simulations  of  basal  and  rhombohedral  twinning  in  a  sap¬ 
phire  single  crystal  with  a  pre-existing  halfpenny-shaped 
notch  are  analyzed.  Conclusions  follow  in  Section  6. 
Regarding  notation,  vectors  and  second-order  tensors  are 
written  in  bold  italic;  scalars  and  components  are  written 
in  plain  italic,  with  summation  applied  over  repeated  indi- 
cial  subscripts. 

2.  Theory 

Only  essential  details  of  the  phase-field  theory  are 
reported  here;  complete  descriptions  are  given  elsewhere 
[9,10].  Let  x  and  X  denote  sufficiently  smooth  spatial  and 
material  coordinates  of  a  body  of  reference  volume  Q, 
related  by  the  differentiable  mapping  x  =  /(A,  t )  that  is 
one-to-one  and  invertible  at  any  fixed  t.  Let  rj(X ,  t)  denote 
the  order  parameter  that  distinguishes  between  the  original 
(parent)  crystal,  the  twin,  and  the  interfacial  boundary 
regions  between  parent  and  twin: 

rj  =  OVA  G  parent,  rj  =  1VA  G  twin, 

0  <  rj  <  1VA  G  boundary.  (2.1) 

The  deformation  gradient  is 

F  =  X%  =  FeF\  (2.2) 

where  V  denotes  the  material  gradient,  FE  accounts  for 
elastic  stretch  and  rotation,  and 

Fn  =  1  +  [<p(rj)y0]s  0  m  (2.3) 

accounts  for  twinning  shear.  Orthogonal  unit  vectors  (in 
material  coordinates)  in  the  twinning  direction  and  normal 
to  twinning  plane  are  s  and  m ,  the  magnitude  of  maximum 
twinning  shear  is  y0?  and  (p(o )  is  an  interpolator 
satisfying  cp( 0)  =  0,  cp(\)  =  1,  cp'( 0)  =  cp'(  1)  =  0,  where 
(■)'  =  d(-)  /  drj .  Defining  CE  =  ( FE)TFE ,  the  local  ratio  of  de¬ 
formed  to  initial  volume  is  J  =  detF=  (detC^)1/2. 

The  total  energy  functional  for  the  body  is 

ti)=  [  [W{F,  rj)  +  /fa,  Vi/)] dG.  (2.4) 

Jq 

The  strain  energy  W  and  interfacial  energy  /  per  unit  refer¬ 
ence  volume  are  of  the  form 
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W=  W[CE (F,ri),  ri],  f  =  f0(r,)  +  k 

:  (Vfy  (8)  Vi/).  (2.5) 

Let  denote  the  boundary  of  Q.  Imposing  the  variational 
principle 

5Y=  <£  {t-5z  +  hdo)dS,  (2.6) 

JdQ 

the  following  local  equilibrium  equations  and  boundary 
conditions  are  derived  [9]: 

V  •  dW/dF  =  V  •  P  =  0, 

f'  -  2V  •  (fcV/7)  +  =  0;  (2.7) 

t  —  F  n ,  h  =  2k  :  (S7ri  <g>  n).  (2.8) 

Traction  per  unit  reference  area  is  /;  conjugate  force  to  the 
order  parameter  is  /z;  outward  normal  to  is  n.  Phase 
equilibrium  in  (2.7)  can  be  written  for  homogeneous  k  as 

/o  =  2k  :  VVf/  +  Ty0<p'- 

t  =  S  :  [C£(s  <8>  m)F"  “*],  S  =  2 dW/dCE.  (2.9) 

In  the  present  work,  attention  is  restricted  to  a  double¬ 
well  potential  /0  =  Arj2(\  —  r\)2,  with  A  a  constant  related 
to  the  barrier  height.  When  isotropic  surface  energy  is 
imposed,  for  which  k  =  /cl,  with  k  a  constant,  then 

/  =  Arj2(  1  -  rj)2  +  k|V/?|2,  k  =  3/7/4, 

=  12T//,  (2.10) 

where  F  and  /  are  equilibrium  energy  per  unit  area  and 
thickness  of  an  unstressed  interface  [9].  Anisotropic  repre¬ 
sentations  of  surface  energy  are  also  considered.  For  exam¬ 
ple,  in  a  material  coordinate  frame  {ej}  with  e\\\s  and  e2\\m , 
then  K\\>  k22  accounts  for  the  increase  in  energy  at  an 
incoherent  twin  boundary  [20]  due  to  the  local  (e.g.  core) 
energy  of  twinning  dislocations  [21],  and  Ki2  =  zc2i  =  0  in 
this  coordinate  system. 

Several  different  elastic  models  are  considered  [10].  For 


compressible  neo-Hookean  behavior 

w  =  tA[(lnJ)2  +  trC£-3]  -  (ilnJ,  (2.11) 

where  2  and  are  usual  isotropic  elastic  constants.  The 
First  Piola-Kirchhoff  stress  is 

P  =  dW/dF  =  FESFn~T 

=  FE[fil  +  (/Un  J  —  n)CE~x]Fn~T .  (2.12) 

For  linear  elastic  behavior, 

W=W[eE(Vu,r,U)  =  l-CUKL(rtK4L,  (2.13) 

where  the  following  geometric  relationships  apply: 

F=l  +  V»~l+J?£+0";  (2.14) 

Fe~1  +  Pe,  eE  =  ^]fiE  +  (pE)J],  CE ~  1  +  2e£;  F"  =  1  +  p";  (2.15) 
P  =  dW/dVu,  i  =  P :  (s 0 m) .  (2.16) 


For  anisotropic  elasticity,  second-order  coefficents  C iJKL 
are  interpolated  between  parent  and  twinned  domains 
using  cp  [9,10].  For  isotropic  elasticity,  Cijkl  =  ^ij^kl  + 
SilSjk),  and  W  and  CIJKL  do  not  depend 
explicitly  on  r\.  The  elastic  driving  force  for  twinning, 
t,  is  the  resolved  shear  stress  on  the  twinning  plane  in  the 
direction  of  twinning  shear. 

Two  different  interpolation  functions  are  also 
considered: 

(p  =  (3  -  2rj)rj2  (polynomial),  (2.17) 

cp  =  1/[1  +  e-2^-a5)]  (exponential).  (2.18) 

Polynomial  function  (2.17)  has  been  used  frequently 
[9,10,22]  and  yields  a  gradual  change  in  cp  over  0  <  rj  <  1; 
the  Fermi-Dirac  function  (2.18),  here  with  k—  15,  provides 
a  steeper  increase  in  cp  at  rj  «  0.5,  as  is  clear  from  Fig.  1. 

Solutions  to  governing  equations  are  obtained  numeri¬ 
cally  using  the  finite-element  method.  The  solution  proce¬ 
dure  [9]  involves  minimization  of  free  energy  functional 
T',  subject  to  possible  boundary  conditions/constraints, 
over  domain  Q,  yielding  the  equilibrium  fields  (/,  //).  Kinet¬ 
ics  and  dissipation  are  not  quantified  explicitly. 

3.  Twinning  under  mode  I  and  II  loading  in  isotropic  solids 

3.1.  Dimensionless  parameters 

Consider  the  phase-field  theory  of  Section  2  applied  to 
an  isotropic  solid.  Dividing  by  the  shear  modulus,  a  nor¬ 
malized  free  energy  functional  becomes 

[  ( W+f)dQ ,  (3.1) 

Jq 

W  =  i  V^(trV«)2  +  (Vi*  —  yo<ps0jw)s :  (Vw-y0<ps(g)m)s,  (3.2) 

f=Arj2(l-r])2  +  Kl2\Vrj\\  (3.3) 

v  =  2/(22  +  2/0,  A  =A/ix  =  12  r/(jd),  k  =  k/^12)  =  3/7(4/*/).  (3.4) 

Notation  (*)s  denotes  the  symmetric  part  of  a  second-order 
tensor.  Since  ju  >  0,  a  solution  («,  rj)  f°r  a  given  set  of 
boundary  conditions  on  dO  that  minimizes  ¥  also  mini¬ 
mizes  P;  this  could  be  a  local  (metastable)  or  global  (sta¬ 
ble)  minimum  energy  configuration. 


Fig.  1.  Phase-field  interfacial  interopolation  functions  (2.17)  and  (2.18). 
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Let  r0  denote  a  constant  reference  value  of  twin  bound¬ 
ary  surface  energy  per  unit  area,  and  let  A0  and  k0  denote 
the  corresponding  values  of  parameters  in  (3.4)  for  fixed 
H  =  Ho  and  fixed  /.  Then 

r/r0  =  A/A0  =  k/k0.  (3.5) 

In  all  subsequent  analysis  /  is  held  fixed.  In  what  follows, 
representative  values  of  F0  =  100  mJ  m-2  and 
Ho  =  25  GPa,  with  /  =  1.0  nm  [9,10,22],  are  used  to  establish 
A o  and  kq.  During  simulations,  the  effect  of  surface  energy  is 
then  studied  by  varying  r/r0.  In  most  simulations,  isotropic 
surface  energy  is  assumed,  but  in  some  (2-D)  simulations, 
Kn  =  (xk22  =  ock  is  prescribed,  where  a>l  accounts  for 
incoherent  boundary  energy  as  noted  in  Section  2.  (For  a 
twin  system  oriented  differently  from  s\\e\  and  m\\e2,  k  is  ro¬ 
tated  as  a  second-order  tensor.)  Accordingly,  for  a  given  set 
of  boundary  conditions  and  prescribed  initial  crystal  orien¬ 
tation  (5,  m ),  solutions  thus  obtained  depend  on  the  param¬ 
eter  set  (y0,  v,  r/r0,  oc)  and  choice  of  interpolation  function 
cp.  The  same  arguments  apply  for  the  isotropic  neo-Hook- 
ean  elastic  model  of  Section  2. 

Plane-strain  simulations  of  an  initially  square  domain  Q 
of  size  2 a  x  2 a  are  reported  in  Section  4.  The  domain  con¬ 
tains  a  pre-existing  straight  edge  crack  of  length  a  and 
thickness  2r0,  with  a  rounded  tip  of  radius  r0.  The  crack 
is  assigned  a  finite  radius  for  two  reasons:  (i)  a  perfect  slit 
crack  would  result  in  truly  singular  stress  fields  at  the  tip 
that  cannot  be  fully  resolved  with  conventional  finite  ele¬ 
ments;  and  (ii)  in  a  real  material  a  finite  r0  on  the  order 
of  or  exceeding  a  lattice  parameter  is  expected  since  the 
separation  distance  across  opposite  faces  of  the  crack 
should  exceed  a  cut-off  distance  for  interatomic  cohesive 
forces  that  would  otherwise  result  in  traction  across  oppos¬ 
ing  faces  of  the  crack. 

Let  the  origin  of  reference  coordinate  systems  ((r,  0)  in 
polar  form  or  ( X ,  Y)  in  rectangular  form)  be  located  at 
the  crack  tip.  Boundary  conditions  are  imposed  as  follows. 
The  crack  surface  E  (6  =  ±n  rad,  0  <  r  <  a)  is  free  of  trac¬ 
tion  (t  =  0).  Neumann  conditions  h  =  0  for  conjugate  force 
to  the  order  parameter  are  assigned  along  all  of  dQ.  Along 
external  boundary  dQ\E  corresponding  to  X,  Y  =  ±a,  dis¬ 
placements  «(r,  6)  corresponding  to  the  K  field  for  pure 
mode  I  or  mode  II  loading  [23]  are  imposed.  For  mode  I: 


ux  =  2A [ar/ (2n)\ 1/2  cos(0/2)  [1  -  2v  +  sin2(0/2)],  (3.6) 

uy  =  2A[ar/(2n)\111  sin(0/2)[2  —  2v  —  cos2(0/2)],  (3.7) 

A  =Kl/{2nal/2).  (3.8) 

For  mode  II: 

iix  =  2 A[ar/ (27i)]^2  sin(0/2)[2  —  2v  +  cos2(0/2)],  (3.9) 

iiY  =  —  2A  [ar / (In)]1^2  cos(0/2)[l  —  2v  —  sin2 (0/2)] ,  (3.10) 
A  =Kll/{lHal/2).  (3.11) 


For  both  modes,  the  orientation  of  the  twin  system  (s,  m)  is 
such  that  the  resolved  shear  stress  t  of  (2.16)  is  maximum 


according  to  the  linear  elastic  solution  [23].  For  mode  II, 
s  is  simply  oriented  in  the  sense  of  positive  r  along  6  =  0. 
For  mode  I,  s  is  oriented  in  the  sense  of  positive  r  along 
6  =  1.22  rad. 

During  phase-field  simulations,  A  is  increased  incremen¬ 
tally.  For  each  increment,  the  domain  is  seeded  with  a  small 
twin  nucleus  at  r  <  r0.  Displacement  boundary  conditions 
are  updated  according  to  the  mode  of  loading  via  (3.6), 
(3.7),  (3.8)  or  (3.9),  (3.10),  (3.11),  and  then  the  equilibrium 
solution  («,  t])  in  Q  is  obtained  through  energy  minimization 
using  the  finite-element  method.  If  the  driving  force  for 
twinning  is  insufficient,  then  the  initial  nucleus  will  disap¬ 
pear,  and  the  equilibrium  solution  includes  t]  =  OVA  e  Q. 
Otherwise,  at  a  threshold  load  parameter  A  =  AT,  a  stable 
twin  will  appear  at  the  crack  tip  (r  =  0).  With  further 
increasing  A  >  AT,  the  twin  will  grow  in  length  and/or  thick¬ 
ness  until  it  interacts  with  the  external  boundary  dQ\E. 

According  to  linear  elastic  fracture  mechanics,  crack 
extension  (i.e.  cleavage)  will  occur  if  the  applied  stress 
intensity  factor  or  corresponding  strain  energy  release  rate 
exceeds  a  threshold  for  a  particular  material  and  loading 
mode: 

Kyu  ^  Kc  Gy  ii  ^  Gc  =>  fracture,  (3-12) 

where  the  fracture  surface  energies  are 
Gl/u=K21/u(l-v)/(2^, 

Gc  =  K2C(\  —  v)/ (2/x),  (3.13) 

and  here  no  distinction  has  been  made  in  notation  among 
threshold  fracture  energies  Gc  for  different  modes.  For 
comparison,  dimensionless  twinning  and  fracture  parame¬ 
ters  associated  with  the  normalized  strain  energy  required 
for  either  mechanism  can  be  constructed: 

rT  =  rT/(nl)  »  a(\  -  v)A2T/l,  Gc  =  Gc/(nl).  (3.14) 

The  following  criteria  then  emerge  that  predict  either  crack 
extension  or  twin  emission  from  the  crack  tip: 

2Fy/Gc  =  2rT/Gc  »  1  =>  fracture, 

2rT/Gc  =  2 rT/Gc  <C  1  =>  twinning.  (3.15) 

Dimensionless  FT  can  be  interpreted  as  an  inverse  measure 
of  the  “twinnability”  of  a  given  material  subjected  to  mode 
I  or  mode  II  loading,  with  smaller  values  of  FT  denoting  an 
increased  tendency  for  crack  tip  twinning.  The  factor  of 
two  arises  because,  in  the  usual  convention  of  fracture 
mechanics,  the  strain  energy  release  rate  Gc  is  twice  the 
fracture  surface  energy  Fc.  (In  this  paper,  notation  “F” 
is  associated  generically  with  surface  energy,  66  G”  with 
strain  energy  release.)  When  2 FT  «  Gc,  strain  energy  re¬ 
leased  by  twinning  and  crack  extension  are  comparable, 
and  either  mechanism  could  be  expected  to  occur.  Note 
that  because  A  oc  a~l/2,  imposed  displacements  u  and  twin¬ 
ning  resistance  FT  do  not  depend  on  a ,  which  serves  merely 
as  a  normalization  constant  to  render  these  quantities 
dimensionless. 
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In  phase-field  simulations,  meshes  are  used  with  signifi¬ 
cant  refinement  (element  size  <C/)  near  the  crack  tip  and 
along  the  anticipated  path  of  twin  extension,  such  that 
solutions  are  independent  of  mesh  resolution.  Twin  nucle- 
ation  depends  strongly  on  local  fields  near  r  —  0  but  not 
strongly  on  a\  further  increasing  the  domain  size  above 
all  =  50  did  not  significantly  affect  initial  stages  of  twin¬ 
ning,  but  as  the  twin  grows  and  approaches  the  boundary, 
the  choice  of  a  necessarily  affects  the  solution.  Solutions 
can  be  modestly  dependent  on  notch  radius  for  small  r0, 
so  two  choices  of  r0  are  explored.  Tables  1  and  2  list  param¬ 
eters  investigated  in  simulations  reported  subsequently  in 
Section  4.  The  typical  (i.e.  usual)  parameter  set  referred 
to  as  “linear  elastic”  in  subsequent  figures  is  given  in 
Table  1.  Normalized  twin  boundary  energy  is 
T0  =  r/(nl).  Deviations  from  these  parameters  referred 
to  in  some  later  figures  are  explained  in  Table  2.  In  partic¬ 
ular,  the  value  of  a  =  100  follows  from  Refs.  [20,21]. 

3.2.  Numerical  results 

Figs.  2  and  3  show  characteristic  results  for  mode  I  and 
mode  II  loading,  respectively.  The  undeformed  finite-ele¬ 
ment  (FE)  mesh  is  shown  in  part  (a)  of  each  figure.  Twin 
morphology  and  a  stress  component — tensile  stress  for 
mode  I,  shear  stress  for  mode  II — are  shown  in  parts  (b) 
and  (c),  corresponding  to  a  load  increment  exceeding 
nucleation,  i.e.  A  >  AT.  In  each  case,  the  twin  (i.e.  stress- 
free  shear  eigenstrain  yd 2)  relieves  much  of  the  stress  that 
would  otherwise  be  large  as  r  — >  0  in  an  elastic  medium 
without  a  twin.  Twin  growth  to  the  boundary  dQ  is  inhib¬ 
ited  by  the  displacement  boundary  conditions. 

Twin  morphologies  for  various  simulations  of  mode  I 
deformation  are  compared  at  the  same  load  increment 
A  =  0.04  >  At,  i.e.  at  the  same  imposed  Kh  in  Fig.  4. 
Because  differences  in  twin  size  and  shape  are  small,  it  is 
concluded  that  model  predictions  of  fully  formed  twin 
morphology  are  not  sensitive  to  choice  of  linear  or  nonlin¬ 
ear  (i.e.  neo-Hookean)  elastic  model,  choice  of  interpola¬ 
tion  function  (2.17)  or  (2.18),  twin  boundary  surface 
energy  anisotropy  a,  or  pre-existing  thickness  of  the  crack 
or  notch.  However,  as  discussed  later,  twin  nucleation  is 
affected  by  these  choices  in  some  cases. 

Fig.  5  shows  the  effects  of  dimensionless  material  prop¬ 
erties  on  twinning  resistance  2FT  of  (3.14)  under  mode  I 
loading.  Recall  from  (3.15)  that  this  resistance  can  be  com¬ 
pared  with  Gc  to  predict  whether  twinning  or  cleavage 
crack  extension  would  be  energetically  favorable,  with 
smaller  FT  suggesting  a  greater  tendency  for  twinning  at 
the  crack  tip.  Each  data  point  on  each  piecewise  linear 
curve  in  Fig.  5  represents  the  result  of  a  different  phase-field 

Table  1 

Basic  simulation  parameters. 

Descriptor  Elasticity  model  (p  r0/l  a  To 

Linear  elastic  Linear  isotropic  Polynomial  2  1  4  x  10-3 


Table  2 

Other  simulation  parameters. 

Descriptor 

Difference  from  basic  parameters 

Neo-Hookean 

Neo-Hookean  elastic  energy 

Exponential  interpolant 

Exponential  cp 

Anisotropic  k 

a  =  kii/k22  =  100 

Thin  notch 

r0/l  =0.5 

simulation  in  which  A  (i.e.  K{)  is  increased  in  increments 
of  10-3  from  A  =  0  to  the  condition  for  which  a  twin  or 
twin  nucleus  is  first  observed  at  A  =  AT. 

Effects  of  twinning  shear  y0  on  crack  tip  twinnability  are 
shown  in  Fig.  5a,  where  discrete  values  of  y0  =  (0.1,  0.3, 
0.5,  0.7,  1.0)  have  been  prescribed  in  simulations  incorpo¬ 
rating  linear  elastic  or  neo-Hookean  strain  energy  density 
W.  Twinning  shear  significantly  affects  nucleation.  A  min¬ 
imum  of  twinning  resistance  FT  is  predicted  at  y0  =  0.3  for 
each  model.  For  y0  <  0.3,  resistance  increases  since  the 
eigenstrain  does  not  reduce  elastic  energy  so  much.  Nucle¬ 
ation  resistance  also  increases  for  y0  >  0-3,  presumably 
because  the  applied  Kx  field  must  be  sufficiently  strong  such 
that  a  large  eigenstrain  relieves  the  elastic  stress  field 
induced  by  the  crack. 

The  effects  of  Poisson’s  ratio  v  on  twin  nucleation  resis¬ 
tance  are  comparatively  small,  as  shown  in  Fig.  5b  for  val¬ 
ues  of  v  =  (0.05,  0.15,  0.25,  0.35,  0.45).  The  low  influence  of 
Poisson’s  ratio  on  twinning  found  here  agrees  with  conclu¬ 
sions  of  a  previous  linear  elastic  analysis  [24].  Neo-Hook¬ 
ean  elasticity  is  more  sensitive  than  linear  elasticity  to  v, 
as  expected  considering  the  nonlinear  compressibility 
inherent  in  (2.11). 

As  shown  in  Fig.  5c,  twin  boundary  energy  F0  strongly 
affects  twinnability,  with  resistance  FT  increasing  with 
increasing  F/r0  in  all  cases.  Discrete  values  F/r0  =  (0.5, 
0.75,  1,  1.5,  2)  have  been  probed.  Twinning  resistance 
increases  relative  to  the  linear  elastic  case  when  the  expo¬ 
nential  interpolator  of  (2.18),  anisotropic  surface  energy 
(a  =  100)  or  a  thinner  notch/crack  is  used.  Differences 
increase  as  the  ratio  F/r0  increases. 

Fig.  6  shows  effects  of  dimensionless  material  properties 
on  twinning  resistance  2FT  of  (3.14)  under  mode  II  loading, 
and  is  analogous  to  results  for  mode  I  of  Fig.  5.  Again, 
each  data  point  represents  the  result  of  a  different  phase- 
field  simulation  in  which  A  (here  proportional  to  Ku)  is 
increased  in  increments  of  1  x  10-3  from  A  =  0  to  the  con¬ 
dition  for  which  a  twin  or  twin  nucleus  is  first  observed  at 
A  =  Ay. 

The  effects  of  twinning  shear  y0  on  crack  tip  twinnability 
are  shown  in  Fig.  6a,  where  discrete  values  of  y0  =  (0.1,  0.3, 
0.5,  0.7,  1.0)  have  been  prescribed  in  simulations  incorpo¬ 
rating  linear  elastic  or  neo-Hookean  strain  energy  W. 
Twinning  shear  significantly  affects  mode  II  nucleation, 
as  was  the  case  for  mode  I.  Here,  a  minimum  of  twinning 
resistance  FT  is  predicted  at  y0  =  0.5  for  each  model  in 
mode  II  loading,  which  exceeds  the  minimum  associated 
y0  =  0.3  observed  for  mode  I. 
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(a)  (b)  (c) 


Fig.  2.  Mode  I  loading  of  isotropic  elastic  solid  [y0  =  \ ,  v  =  J ,  T  =  F0] :  (a)  finite-element  mesh;  (b)  order  parameter  at  A  =  0.05;  (c)  tensile  normal  stress 
(P]  i  =  PYy )  at  A  =  0.05.  The  origin  of  the  (X,  Y)  coordinate  system  is  at  the  crack  tip,  with  positive  X  downward  and  positive  Y  to  the  right.  For  polar  (r, 
0)  coordinates,  0  is  measured  counterclockwise  from  the  positive  X-axis. 


(a)  (b)  (c) 


Fig.  3.  Mode  II  loading  of  isotropic  elastic  solid  [y0  =  \ ,  v  =  \ ,  T  =  F0] :  (a)  finite-element  mesh;  (b)  order  parameter  at  A  =0.03;  (c)  shear  stress 
(P 12  =  Pxy=  Pyx)  at  A  =  0.03.  The  origin  of  the  (X,  Y)  coordinate  system  is  at  the  crack  tip,  with  positive  X  downward  and  positive  Y  to  the  right.  For 
polar  (r,  0)  coordinates,  6  is  measured  counterclockwise  from  the  positive  X-axis. 


(a)  (b)  (c) 


(d)  (e) 


Fig.  4.  Order  parameter  (twin  morphology)  for  mode  I  loading  of  isotropic  elastic  solid  at  A  =  0.04  [y0  =  \ ,  v  =  \ ,  P  =  F0] :  (a)  linear  elastic;  (b)  neo- 
Hookean;  (c)  exponential  interpolant;  (d)  anisotropic  k;  (e)  thin  notch. 

The  effects  of  Poisson’s  ratio  v  on  twin  nucleation  resis¬ 
tance  under  mode  II  loading  are  comparatively  small,  as 
shown  in  Fig.  6b.  Twinning  resistance  for  neo-Hookean 
elasticity  is  again  more  sensitive  than  linear  elasticity  to  v. 


As  shown  in  Fig.  6c,  twin  boundary  energy  F0  strongly 
affects  twinnability,  with  resistance  FT  increasing  with 
increasing  F/r0  in  all  cases,  as  was  observed  for  mode  I. 
Twinning  resistance  increases  relative  to  the  linear  elastic 
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r/r0 

(c) 


Fig.  5.  Normalized  twin  nucleation  energy  FT  under  mode  I  loading  for  (a)  variable  twinning  shear  y0;  (b)  variable  Poisson’s  ratio  v;  and  (c)  variable  twin 
boundary  surface  energy  r.  A  low  value  of  FT  correlates  with  a  low  value  of  applied  Ki  needed  to  initiate  a  twin  at  the  crack  tip. 


case  when  the  exponential  interpolator  of  (2.18)  or  aniso¬ 
tropic  surface  energy  (a  =  100)  is  used,  though  the  latter 
influences  the  results  modestly  and  only  for  r/r0  >  1.  Con¬ 
trary  to  results  for  mode  I  loading,  a  thinner  notch  reduces 
rather  than  increases  twin  resistance  under  mode  II.  Differ¬ 
ences  among  cases  in  Fig.  6c  increase  as  F/r0  increases. 

Comparing  results  in  Figs.  5  and  6,  numerical  values  of 
resistance  to  twinning  FT  under  mode  I  loading  tend  to 
exceed  those  under  mode  II  by  a  factor  of  the  order  of  2. 
For  example,  for  linear  elasticity  with  the  parameter  set 
[y0  =  y ,  v  =  f ,  T  =  r0],  phase-field  simulations  predict 
rT  «  0.007  for  mode  II  and  FT  «  0.011  for  mode  I.  This 
result  is  not  unexpected  since  a  minimum  twinning  resis¬ 
tance  would  be  associated  with  the  geometry  of  mode  II 
loading  in  Fig.  3,  i.e.  with  the  twin  system  aligned  perfectly 
with  maximum  shear  stress  at  the  tip  of  a  mode  II  crack. 

4.  Twinning  under  mode  I  loading  in  single  crystals 

Twinning  under  mode  I  loading  of  a  square  domain  with 
pre-existing  edge  crack  is  investigated  next  for  single  crys¬ 
tals  with  relevant  physical  properties.  Plane-strain  simula¬ 
tions  on  a  sample  of  dimensions  identical  to  that  of 
Section  3  are  reported.  In  this  2-D  idealization,  only  one 
twin  system  is  permitted  to  be  active  in  any  simulation, 
and  the  crack  propagation  direction  in  a  particular 
cleavage  plane  is  chosen  such  that  crack  opening  is  in  the 
plane  0  =  ±7i  rad,  i.e.  the  pre-existing  crack  is  along 


Y  =  0,X  <  0.  Details  regarding  properties  and  results  are 
reported  in  Table  3;  corresponding  discussion  for  each 
material  follows  next.  As  will  become  clear  later,  “Model” 
in  Table  3  designates  the  type  of  elasticity  model  and/or 
twin  boundary  surface  energy  representation,  with  “isotro¬ 
pic”  referring  to  isotropic  linear  elasticity  and  isotropic  sur¬ 
face  energy  (a=l),  “aniso.  IF”  referring  to  anisotropic 
linear  elasticity  and  isotropic  surface  energy,  and 
“a  =  100”  denoting  isotropic  linear  elasticity  with  aniso¬ 
tropic  surface  energy. 

4.1.  Calcite 

Calcite  is  a  soft  mineral  of  trigonal  symmetry  whose 
pure  crystals  are  transparent.  Calcite  twins  readily,  with  lit¬ 
tle  or  no  plastic  slip,  under  concentrated  surface  loading. 
The  preferred  twin  system  is  e+,  with  relatively  large  shear 
yo  =  0.694  and  geometry  (100) {011}  in  rhombohedral 
pseudocell  notation  [25].  Calcite  also  cleaves  easily  on  the 
natural  rhombohedral  planes  (i.e.  {100}  planes)  of  its 
primitive  unit  cell,  equivalent  to  {1  Ol  1}  planes  in  the  hex¬ 
agonal  notation  of  Refs.  [26,27].  In  the  present  simulations, 
a  2-D  projection  is  required,  where  0  =  0.89  rad  is  the 
resulting  orientation  of  the  e+  twin  system  that  receives 
the  maximum  stress  t  under  mode  I  loading  of  a  cleavage 
plane.  Cleavage  surface  energy  entering  Gc  in  Table  3  is 
obtained  from  experiments  [27].  Properties  associated  with 
twinning  and  elasticity  are  from  Ref.  [10]  and  references 
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Yo 

(a) 


(b) 


r/r0 

(c) 


Fig.  6.  Normalized  twin  nucleation  energy  rT  under  mode  II  loading  for  (a)  variable  twinning  shear  y0;  (b)  variable  Poisson’s  ratio  v;  (c)  variable  twin 
boundary  surface  energy  f.  A  low  value  of  FT  correlates  with  a  low  value  of  applied  Ku  needed  to  initiate  a  twin  at  the  crack  tip. 


Table  3 

Single-crystal  properties  and  results  of  phase-field  simulations. 


Material 

V 

p  (GPa) 

Crack 

Gc 

Twin  system 

7o 

2r0 

Model 

2/VGc 

Prediction 

Calcite 

0.30 

40 

(1011) 

0.017 

e+ 

0.694 

0.0091 

aniso.  W 

2.53 

Fracture 

Isotropic 

2.96 

Fracture 

Sapphire 

0.23 

167 

(0001) 

0.48 

Rhomb.  (R) 

0.202 

0.0015 

Isotropic 

0.02 

Twinning 

Max.  t 

0.1 

0 

0 

II 

0.11 

Twinning 

(1012) 

0.072 

Basal  (B) 

0.635 

0.0089 

Isotropic 

0.62 

Either 

Max.  t 

0.1 

0 

0 

II 

0.65 

Either 

Mg 

0.28 

19 

(0001) 

382 

[1  Ol  1]  (1 0 1 2) 

0.130 

0.0121 

Isotropic 

2  x  10“4 

Twinning 

therein.  Both  isotropic  and  anisotropic  elastic  models  are 
investigated  (i.e.  forms  of  W).  For  the  former,  the  Voigt 
elastic  constants  shown  in  Table  3  apply.  For  the  latter, 
values  (Cn  =  165,  Ci2  =  65,  Ci3  =  62,  Ci4  =  —23, 
C33  =  89,  C44  =  37  GPa)  from  Ref.  [10]  are  used.  For  either 
elastic  model,  the  displacement  field  of  (3.6)  and  (3.7)  is 
applied  by  incrementally  increasing  A  or,  equivalently, 
Ki,  this  is  an  approximation  of  the  true  Kx  field  when  aniso¬ 
tropic  elasticity  is  used. 

Results  in  Fig.  7a  and  b  show  the  twin  at  loading  soon 
after  nucleation,  i.e.  at  A  =  0.044  >  AT.  Nucleation 
occurs  first  for  the  anisotropic  model,  but  the  orientation 
(0  «  0.9  rad)  and  shape  of  the  twin  are  similar  in  each  case. 
A  secondary  twin  belonging  to  the  same  twin  system  begins 
to  form  at  a  larger  applied  K\  field,  as  shown  in  Fig.  7c.  For 


anisotropic  and  isotropic  models,  2.5  <  2 rT/Gc  <  3, 
meaning  that  crack  extension  is  preferred  over  twinning 
according  to  criterion  (3.15).  This  result  is  in  qualitative 
agreement  with  tensile  fracture  experiments  [26,27]  that 
report  no  evidence  of  twinning.  Parting  fractures  induced 
by  twins  in  calcite  have  also  been  noted  [28].  These  model 
results  do  not  contradict  the  possibility  of  twins  induced  by 
defects  during  other  modes  of  loading,  e.g.  under  spherical 
indentation,  samples  with  visible  surface  cracks  are  known 
to  twin  more  easily  than  those  without  [5]. 

4.2.  Sapphire 

Sapphire,  also  known  as  corundum  or  single-crystal  alu¬ 
mina,  is  a  hard  ceramic/mineral  that,  like  calcite,  is  of  tri- 
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Fig.  7.  Order  parameter  for  mode  I  cleavage  of  calcite  single  crystal:  (a)  anisotropic  elasticity,  A  =  0.044;  (b)  isotropic  elasticity,  A  =  0.044;  (c)  isotropic 
elasticity,  A  =  0.1. 


gonal  elastic  symmetry  and  can  be  transparent.  As 
reviewed  in  Refs.  [4,10],  the  twin  systems  are  rhombohedral 
(R)  with  Miller  indices  <1 10T){1  TO 2}  in  the  structural 
unit  cell,  and  basal  (B)  with  Miller  indices  <1  T 00){000 1} 
in  the  structural  unit  cell.  Twinning  shear  for  R  systems 
(y0  =  0.202)  is  less  than  that  for  B  systems  (y0  =  0.635). 
Surface  energies  for  cleavage  on  rhombohedral,  prismatic 
and  basal  planes  have  been  reported,  with  R  planes  the 
most  likely  to  cleave  and  B  planes  difficult  to  fracture 
[29].  In  various  phase-field  simulations  reported  next,  either 
a  B  twin  system  or  a  R  twin  system  is  active,  with  a  mode  I 
crack  located  on  one  of  several  possible  planes.  Specifically, 
the  four  cases  reported  in  Table  3  correspond  to  (i)  R  twin¬ 
ning  induced  by  a  basal  plane  crack  (s  at  6  =  1.07  rad);  (ii) 
R  twinning  induced  by  a  noncrystallographic  plane  crack 
that  produces  maximum  t  (5  at  0  =  1.22  rad);  (iii)  B  twin¬ 
ning  induced  by  a  rhombohedral  plane  crack  (s  at 
9  =  1.07  rad);  and  (iv)  B  twinning  induced  by  noncrystallo¬ 
graphic  plane  crack  that  produces  maximum  t  (5  at 
0=1.22  rad).  Cleavage  surface  energies  entering  Gc  in 
Table  3  are  obtained  from  [29].  For  cases  (i)  and  (iii),  iso¬ 
tropic  surface  energy  is  used.  For  cases  (ii)  and  (iv),  the 
effects  of  anisotropic  k  associated  with  incoherent  twin 
boundary  energy  are  explored  by  setting  a  =  100  [21,20]. 
In  all  cases,  isotropic  elasticity  is  imposed  with  Voigt  elastic 
constants,  noting  from  previous  work  [10]  that  the  effects  of 
elastic  anisotropy  are  small  in  sapphire;  elastic  anisotropy 
is  also  investigated  later  in  3-D  simulations  in  Section  5, 
confirming  this  assertion. 


Results  for  cases  (i)  and  (iii)  are  shown  in  Fig.  8a  and  b 
at  A  >  AT.  The  basal  twin  (Fig.  8a)  nucleates  at  a  larger  A 
and  is  thinner  than  the  rhombohedral  twin  (Fig.  8b).  Twin¬ 
ning  resistance  is  compared  with  fracture  energy  in  Table  3. 
Since  2FT  <C  Gc  for  cases  (i)  and  (ii)  involving  R  twinning, 
this  twinning  mode  is  preferred  over  mode  I  crack  exten¬ 
sion.  On  the  other  hand,  2 FT  is  smaller  than  Gc,  but  not 
significantly  so,  for  cases  (iii)  and  (iv)  involving  basal  twin¬ 
ning.  It  follows  that  basal  twinning  is  possible  in  such 
cases,  but  crack  extension  is  also  likely,  considering  possi¬ 
ble  sources  of  uncertainty  in  the  phase-field  model/param¬ 
eters  and  local  variations  in  micro  structure  (e.g.  defects  or 
impurities)  inherent  in  real  experimental  samples.  Predic¬ 
tions  are  in  qualitative  agreement  with  experiments.  Specif¬ 
ically,  in  cleavage  experiments  [29],  basal  fracture  was 
found  much  more  difficult  to  induce  than  rhombohedral 
fracture.  In  recovered  specimens  fractured  by  bending 
[30]  on  unidentified  planes,  numerous  thicker  R  twins  were 
found,  and  fewer  thinner  B  twins  were  observed.  The 
thicker  predicted  shape  of  the  R  twin  relative  to  the  B  twin 
is  evident  in  Fig.  8;  it  has  been  noted  elsewhere  [28]  that 
twin  systems  with  larger  y0  are  prone  to  yield  thinner  twins. 
The  presence  of  both  kinds  of  twins  has  been  reported  in 
shock  compression  experiments  on  alumina  poly  crystals  [3]. 

4.3.  Magnesium 

Magnesium  is  a  moderately  ductile  metal  with  hexago¬ 
nal  crystal  structure.  A  number  of  slip  and  twin  systems 


(a)  (b) 


Fig.  8.  Order  parameter  for  mode  I  cleavage  of  sapphire  single  crystal:  (a)  basal  twinning  and  rhombohedral  cleavage,  A  =  0.057;  (b)  rhombohedral 
twinning  and  basal  cleavage,  A  =  0.046. 
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Fig.  9.  Order  parameter  for  mode  I  basal  plane  cleavage  of  magnesium  crystal:  (a)  twin  nucleation,  A  =  0.044;  (b)  tensile  twinning,  A  =0.1. 


have  been  identified;  the  twin  system  investigated  here  is 
the  dominant  inelastic  mechanism  observed  in  single  crys¬ 
tals  stretched  along  [0001]:  the  system  (1  Ol  1 ) {T 0 1 2}, 
with  relatively  small  shear  y0  =  0.1295.  Elastic  anisotropy 
is  very  low  in  Mg;  Voigt  isotropic  elastic  constants  [9]  are 
used,  along  with  the  isotropic  twin  boundary  surface 
energy  listed  in  Ref.  [9].  Although  various  cleavage  modes 
in  single  crystals  have  been  reported  [31],  quantitative  val¬ 
ues  of  fracture  surface  energies  are  not  evident  in  the  exist¬ 
ing  literature;  however,  analysis  suggests  that  the  fracture 
energies  of  prismatic  and  basal  planes  should  be  approxi¬ 
mately  equal  [32].  In  the  present  work,  a  pre-existing  edge 
crack  on  the  basal  plane  is  modeled,  where  the  value  of  Gc 


in  Table  3  is  obtained  from  the  macroscopic  fracture 
toughness  of  Mg  polycrystals  [33].  The  most  favorably  ori¬ 
ented  twin  system  is  oriented  with  s  at  0  =  0.75  rad. 

The  predicted  twin  is  shown  in  Fig.  9a  at  A  «  AT  and  in 
Fig.  9b  at  A  >  Ay.  The  rounded  shape  of  the  twin  nucleus 
in  Fig.  9a  is  in  qualitative  agreement  with  previous  theoret¬ 
ical  studies  [9,34].  The  symmetric  double-twin  morphology 
in  Fig.  9b  is  similar  to  atomic  simulation  results  of  tensile 
twinning  in  a  Mg  single  crystal  with  a  pre-existing  center 
crack  on  the  basal  plane  [35]  (see  their  Fig.  4).  In  the  pres¬ 
ent  simulations,  the  blunt  shape  of  the  twin(s)  correlates 
with  the  low  value  of  y0  in  Mg.  Twinning  resistance  is  com¬ 
pared  with  fracture  energy  in  Table  3.  Since  2Tt  <C  Gc, 


Fig.  10.  Order  parameter  for  direct  shear  loading  (y  =  1)  of  sapphire  single  crystal  with  halfpenny-shaped  edge  notch:  (a)  basal  twin,  anisotropic  elasticity 
and  anisotropic  surface  energy;  (b)  basal  twin,  isotropic  elasticity  and  isotropic  surface  energy;  (c)  rhombohedral  twin,  anisotropic  elasticity  and 
anisotropic  surface  energy;  (d)  rhombohedral  twin,  isotropic  elasticity  and  isotropic  surface  energy. 


J.D.  Clayton,  J.  Knap  /  Acta  Materialia  61  (2013)  5341-5353 


5351 


(a)  (b)  (c)  (d)  (e) 


Fig.  12.  Twin  nucleation  and  growth  along  mid-plane  Y  =  0  for  direct  shear  loading  of  sapphire  single  crystal,  anisotropic  model,  rhombohedral  twin:  (a) 
y  =  0.2;  (b)  y  =  0.4;  (c)  y  =  0.6;  (d)  y  =  0.8;  (e)  y  =  1.0. 


twinning  is  preferred  over  mode  I  crack  extension  for  the 
present  boundary  conditions. 

The  analytical  model  of  Ref.  [8]  suggests  the  following 
criterion  for  twinning  vs.  crack  extension: 

X  ■  C/t//c)1/2  >  1  =>  fracture, 

X '  C/t//c)1/2  <  1  =>  twinning,  (4.1) 

where  the  dimensionless  parameter  %  depends  on  load 
direction,  crystal  structure  (e.g.  c/a  ratio  in  hexagonal  met¬ 
als)  and  anisotropic  elastic  constants,  and  /T  and  fc  are 
energies  associated  with  “total  inelastic  resistance”  against 
twin  and  crack  extension,  respectively.  For  basal  plane 
cleavage  and  tensile  twinning,  a  value  of  %=  1.66  is  re¬ 
ported  for  Mg  [8].  Values  of  /T  and  fc  have  not  been  re¬ 
ported;  the  former  can  be  deduced  from  the  present 
results  if  fc=\  Gc  is  assumed.  Squaring  both  sides  of 
(4.1)  and  comparing  with  (3.15)  yields 

h  =  rT/x2  =  vITt/x2  «  0.26  J/m2,  (4.2) 

which  is  significantly  larger  than  the  twin  boundary  surface 
energy  F  =  0.12  J/m2.  Such  a  result  reinforces  the  notion 
that  total  energetic  resistance  to  twinning  depends  on  other 
factors  besides  F  alone. 

5.  Twinning  in  a  notched  single  crystal:  3-D  simulations 

Three-dimensional  simulations  demonstrate  how  the 
phase-field  model  can  be  applied  to  predict  twinning  under 
complex  stress  states.  In  the  present  simulations  of  single¬ 
crystal  sapphire,  attention  is  restricted  to  a  single  potentially 
active  twin  system.  Boundary  conditions  imposing  direct, 
intense  shear  strain  resolved  on  this  system  (discussed  in 


detail  below)  are  such  that  only  one  system  would  be 
expected  to  be  active.  Various  simulations  consider  either 
rhombohedral  (R)  or  basal  (B)  twinning,  with  isotropic  or 
anisotropic  material  models.  For  the  former,  isotropic  twin 
boundary  surface  energy  is  also  used.  For  the  latter,  the  tri¬ 
gonal  elastic  constants  (Cn  =  500,  C12  =  168,  C13  =  121, 
C14  =  24,  C33  =  502,  C44  =151  GPa)  from  Ref.  [10]  are  used, 
and  anisotropic  twin  boundary  energy  is  imposed  with 
a  =  100.  The  remaining  material  properties  have  already 
been  discussed  in  Section  4.2  in  the  context  of  Table  3. 

Consider  here  is  a  cube  of  material  with  initial  dimen¬ 
sions  4  ax  4  ax  4a ,  where  a/l=  10.  Six  faces  are  labeled 
±X,  ±Y,  =LZ,  where  the  unit  normal  of  each  face  is 
aligned  parallel  to  the  corresponding  axis  in  a  global 
Cartesian  coordinate  system  with  origin  at  the  center  of 
the  cube.  A  half-cylinder  of  radius  a  and  height  21  is 
extracted  from  the  —X  face  of  the  cube  along  the  mid¬ 
plane  Y  —  0.  This  can  be  interpreted  as  a  pre-existing, 
halfpenny-shaped  notch  or  edge  crack.  Displacement 
boundary  conditions  are  applied  to  the  —X  face  and  cre¬ 
ate  a  region  of  intense  shear  deformation  of  magnitude  y 
over  the  region  —  /  <  Y<  /,  similar  to  conditions  explored 
in  Ref.  [36]  for  modeling  slip,  or  to  what  might  be 
observed  in  the  early  stages  of  a  dynamic  Kalthoff  exper¬ 
iment.  Specifically,  face  —X  is  held  fixed  for 
—2a  <  Y  ^  —l  and  displaced  rigidly  in  the  +X  direction 
for  /  ^  Y  ^  2a.  The  opposite  +X  face  is  held  fixed 
{u  =  0),  and  lateral  faces  =b  Y,  =LZ  are  traction  free.  All 
surfaces  comprising  dQ  (six  cube  faces  and  the  crack  sur¬ 
face)  are  assigned  the  free  ( h  =  0)  Neumann  condition  for 
the  order  parameter,  enabling  possible  twin  nucleation  at 
any  of  these  surfaces. 
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Characteristic  results  are  shown  in  Fig.  10  for  an 
imposed  shear  of  unity  (y  =  1).  Specifically,  B  twinning  is 
depicted  in  Fig.  10a  and  b,  R  twinning  in  Fig.  10c  and  d. 
For  each  mode  of  twinning,  predictions  of  linear  isotropic 
and  anisotropic  elasticity  are  similar  for  the  order  parame¬ 
ter  ( rj )  profile  in  the  region  of  intense  shear.  Anisotropic 
surface  energy  suppresses  formation  of  the  partial  second¬ 
ary  twins  that  emerge  along  the  upper  edge  of  the  —X  face 
in  each  of  the  isotropic  simulations  (Fig.  10b  and  d).  Sim¬ 
ulations  with  neo-Hookean  elasticity  were  also  performed; 
the  results  were  very  similar  to  those  shown  for  linear  iso¬ 
tropic  elasticity  and  are  not  shown.  Consistent  with  the 
results  of  2-D  simulations  in  Section  4.2  and  experimental 
observations  [30],  basal  twinning  is  more  difficult  to  enact 
than  rhombohedral  twinning  under  the  present  direct  shear 
boundary  conditions,  and  B  twins  tend  to  be  thinner  than 
R  twins. 

Order  parameter  contours  along  mid-plane  Y  =  0  are 
shown  in  Figs.  11  and  12,  respectively,  for  B  and  R  twin¬ 
ning  for  incrementally  increasing  applied  deformation  y. 
Some  asymmetry  of  the  twin  boundary  front  is  evident, 
particularly  for  the  B  twin  in  Fig.  11,  a  consequence  of 
anisotropy.  In  these  simulations,  the  semicircular  edge 
crack  does  not  promote  or  inhibit  twinning;  however,  dif¬ 
ferent  boundary  conditions  explored  elsewhere  in  atomic 
simulations  of  shock  compression  [37]  have  demonstrated 
the  possibility  of  R  twinning  induced  at  pre-existing  planar 
cracks  in  sapphire. 

6.  Conclusions 

Twin  emission  from  a  crack  tip  has  been  studied  using 
phase-field  simulations.  A  parameter  associated  with  resis¬ 
tance  to  twin  nucleation  under  mode  I/II  loading  has  been 
derived.  This  parameter  can  be  compared  with  the  fracture 
energy  of  the  material  to  suggest  whether  an  existing  crack 
should  extend  or  a  deformation  twin  should  emerge  and 
grow.  Effects  of  material  properties  and  phase-field  model 
features  on  twinning  resistance  have  been  studied  paramet¬ 
rically,  with  Poisson’s  ratio  and  elastic  nonlinearity  show¬ 
ing  little  effect.  In  contrast,  resistance  to  crack  tip 
twinning  depends  strongly  on  twin  boundary  surface 
energy  and  twinning  shear.  Plane-strain  simulations  of 
twinning  induced  by  a  pre-existing  crack  on  relevant  cleav¬ 
age  planes  in  calcite,  sapphire  and  magnesium  single  crys¬ 
tals  have  been  conducted.  Results  suggest  that  calcite 
should  cleave,  magnesium  should  twin,  and  that  rhombo¬ 
hedral  twinning  is  preferred  to  basal  twinning  in  sapphire, 
all  in  agreement  with  experiment.  Three-dimensional  simu¬ 
lations  of  shear  loading  of  sapphire  demonstrate  a  prefer¬ 
ence  for  rhombohedral  over  basal  twins,  with  the  former 
thicker  in  shape,  in  agreement  with  experiment. 
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